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The objective of this work is to compare a high-order solver with a low-order solver for 
performing large-eddy simulations (LES) of a compressible mixing layer. The high-order method 
is the Wave-Resolving LES (WRLES) solver employing a Dispersion Relation Preserving (DRP) 
scheme. The low-order solver is the Wind-US code, which employs the second-order Roe Physical 
scheme. Both solvers are used to perform LES of the turbulent mixing between two supersonic 
streams at a convective Mach number of 0.46. The high-order and low-order methods are 
evaluated at two different levels of grid resolution. For a fine grid resolution, the low-order 
method produces a very similar solution to the high-order method. At this fine resolution the 
effects of numerical scheme, subgrid scale modeling, and filtering were found to be negligible. 
Both methods predict turbulent stresses that are in reasonable agreement with experimental data. 
However, when the grid resolution is coarsened, the difference between the two solvers becomes 
apparent. The low-order method deviates from experimental results when the resolution is no 
longer adequate. The high-order DRP solution shows minimal grid dependence. The effects of 
subgrid scale modeling and spatial filtering were found to be negligible at both resolutions. For the 
high-order solver on the fine mesh, a parametric study of the spanwise width was conducted to 
determine its effect on solution accuracy. An insufficient spanwise width was found to impose an 
artificial spanwise mode and limit the resolved spanwise modes. We estimate that the spanwise 
depth needs to be 2.5 times larger than the largest coherent structures to capture the largest 
spanwise mode and accurately predict turbulent mixing. 
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speed of sound 

U* 

= 

dimensionless mean velocity 

b 

= 

mixing layer thickness 

Ui 

= 

Favre filtered velocity 

E 

= 

total energy per unit mass 

u 

= 

mean streamwise velocity 
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= 

frequency 

u!u! 
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Reynolds normal stress 

H 
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tunnel height 
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Reynolds shear stress 
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turbulent kinetic energy 
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rms streamwise velocity 
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= 

domain length 

Vrms = 
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rms transverse velocity 
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convective Mach number 

Wrms = yjw'w' 

= 

rms spanwise velocity 

N 

= 

number of grid points 

X i 

= 

Cartesian coordinates 

P 

= 

pressure 

x, y,z 

= 

Cartesian directions 

4 

= 

heat flux 

*ij 

= 

spatially filtered shear stress 

Re 

= 

Reynolds number 

t sgs 

L ij 

= 

subgrid scale shear stress 

T 

= 

temperature 

P 

= 

spatially filtered density 

t 

= 

time 

5* 

= 

displacement thickness 

U 

= 

mean velocity 

Pi 

= 

laminar dynamic viscosity 
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II 

<l 

= 

reference velocity 

Pt 

= 

turbulent dynamic viscosity 

Superscripts 



Subscripts 



~ 


Favre-filtered 

t 


total condition 

- 


Spatially-filtered 

1 


top stream 

SGS 


subgrid scale 

2 


bottom stream 


I. Introduction 

Use of the Reynolds Averaged Navier-Stokes (RANS) equations has become standard practice for aerodynamic 
analysis. However the Reynolds averaging procedure introduces additional correlation terms that must be modeled. 
These models work reasonably well for the equilibrium wall -bounded and free shear layer flows for which they were 
calibrated, but often require retuning of model coefficients for other flows. These models also have severe 
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limitations in predicting unsteady, separated, or transitioning flow. Increasing computational throughput has 
enabled the use of more sophisticated methods such as Direct Numerical Simulations (DNS) and Large-Eddy 
Simulations (LES). These tools are more deterministic and, unlike RANS, are less reliant on tunable constants. 
Consequently, they offer promising physical insights that RANS methods cannot. Despite this, there remains a great 
many issues to ensure the accuracy of LES. 

The primary objective of this work is to compare the benefits and shortcomings of a high-order method and a 
low-order method for performing LES of a compressible mixing layer. The high-order solver is the Wave-Resolving 
LES (WRLES) code [1-3] that employs a Dispersion Relation Preserving (DRP) scheme. The low-order solver is the 
Wind-US code [4] which employs the second-order Roe scheme. It is well known that second-order schemes are 
more convenient to use for complicated geometries due to their small stencil size. Additionally, second-order 
schemes can be easily adapted to unstructured grids whereas high-order schemes frequently are limited to structured 
grids and simple geometries. The advantage of high-order methods is that they can achieve the same level of 
accuracy using significantly fewer grid points compared to the low-order schemes. The work herein sets to examine 
whether low-order solvers, given sufficiently fine enough resolution, can produce the same level of accuracy as 
high-order methods. This effort also attempts to address the current state-of-the-art practice in performing LES of a 
compressible mixing layer. 

Several technologies stand to benefit from an enhanced understanding of compressible mixing in free shear 
flows. Examples include supersonic combustion, exhaust nozzles, and internal flows present in jet engines and 
scramjets, to mention only a few. High speed turbulent flows are often investigated using Reynolds Averaged 
Navier-Stokes (RANS) models, and at times with dilitation-based compressibility corrections. It is widely known, 
however, that these compressibility corrections are empirical at best and do not replicate the actual physics of 
compressible turbulent mixing layers [5-6]. Therefore, the longer term objectives of this work are to better 
understand the fundamentals of compressible turbulent free shear flows using LES and to apply that insight towards 
improving existing RANS turbulence models for engineering applications. 

Brown & Roshko’s [7] classical work on low-speed mixing layers made a key discovery in its noteworthy 
evidence of coherent large-scale spanwise oriented vortical structures. The formation of these large-scale “rollers” is 
triggered by the Kelvin-Helmholtz instability. As the spanwise rollers convect downstream they combine through 
pairing and tearing processes. This work led to a series of experimental studies investigating compressibility effects 
such as Chinzei et al.[8], Papamoschou & Roshko [9], Goebel and Dutton [10], Samimy & Elliot [11], Hall et al. 
[12], and Clemens & Mungal [13]. The experimental studies observed that the most significant effect of 
compressibility on a turbulent mixing layer is suppression of the growth rate relative to an incompressible mixing 
layer at the same velocity and density ratios. The experiments demonstrate that compressibility effects are greatest 
between convective Mach numbers of 0.5 and 1. The particular case examined here is the turbulent mixing between 
two supersonic streams at a moderate convective Mach number of 0.46. Similarly, numerical simulations have been 
performed by several researchers including Foysi and Sarkar [14], Freund et al.[15], and Pantano and Sarkar [16]. 
The numerical simulations showed that dilatation terms are relatively small and not affected by increasing 
convective Mach number, whereas the pressure-strain terms in the Reynolds stress transport decrease monotonically 
with increasing convective Mach number. The redistribution of the Reynolds stresses in turn lowers turbulent 
production, which results in reduced mixing layer growth rate. 

In this work, LES is employed to numerically investigate a reference turbulent mixing layer at moderate 
convective Mach number, 0.46. Two different flow solvers are used. The low-order solver employs a second-order 
upwind scheme while the high-order method implements Bogey & Bailly's [17-18] 11-point DRP scheme. The 
effects of numerical scheme, subgrid scale (SGS) modeling, and filtering are also investigated. Furthermore, a 
previous hybrid RANS -LES work concluded that an insufficient spanwise domain limited the development of the 
full 3D structures [19]. Therefore, the work herein, in addition to its primary purpose of comparing high-order and 
low-order methods, also seeks to analyze the effects of increasing the spanwise extent while maintaining the 
spanwise resolution. This type of parametric study was recently used with the same low-order technique to 
investigate sensitivity of separated bluff body flows to spanwise width [20] . 

Ideally, LES would be performed not only in the mixing region but also along the upstream channels to resolve 
their turbulent boundary layers. Recently, Inoue & Pullin [21] investigated the zero pressure gradient turbulent 
boundary layer through LES over momentum thickness Reynolds numbers of Re d = O(10 3 ) — 0(1O 12 ). This 
required a large number of grid points to resolve the boundary layer. In the present work focusing on the shear layer, 
the Reynolds numbers were Re d = 5.1 x 10 6 and Re Q = 4.7 x 10 6 for the top and bottom channel flows 
respectively. Fully resolving these turbulent boundary layers while also performing LES in the mixing region is 
computationally prohibitive. Therefore, the LES simulations are designed to focus on resolving the shear layer itself. 
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This will also enable a side-by-side comparison between the high-order method WRLES with the low-order method 
Wind-US. 

The paper is organized as follows: A description of the flow is provided in section II. In section III the 
formulation is discussed, while the details of the computational approach are given in section IV. A brief summary 
of the low-order solver Wind-US and the high-order solver WRLES is provided. Details of the grid, boundary 
conditions, and experimental conditions are also presented. Comparisons of the high-order and low-order method 
are presented in section V. The section begins with a general description of the flow and then considers specifics of 
numerical method, SGS modeling, and filtering. The importance of spanwise width is highlighted. The results 
section finishes with a discussion comparing the low-order and high-order method again but at coarser grid 
resolutions. Conclusions are presented in section VI. 

II. Description of the Flow 

The current simulations recreate case 2 of the experimental mixing layer tested by Goebel and Dutton [10]. The 
schematic in Figure 1 shows two supersonic streams separated by a thin splitter plate exiting into a long mixing 
section. The tunnel has a height of 48 mm and a width of 96 mm. The splitter tip separating the two flows has a 
thickness of 0.5 mm (an important length scale affecting vortex shedding) and the entrances of both streams into the 
mixing section have equal heights of 24 mm. The high-speed flow is located on the top with a Mach number of 1.91, 
while the bottom flow is at Mach 1.36. The convective Mach number calculated from the flow conditions is M c = 
0.46, where M c = (U 1 — Ff)/{a 1 + a 2 ). The two flows are pressure matched at 49kPa. The total temperature of the 
top flow is 578 K while the bottom flow is at 295 K. The velocities are U 1 = 700 m/s and U 2 = 399 m/s. 

The experiments used Laser Doppler Velocimetry (LDV) to measure velocities in the mixing section. The total 
viewing length of the optical instrumentation was 500 mm. Transverse velocity profiles at several streamwise 
locations were measured. The Reynolds normal stresses in the streamwise and transverse directions were measured; 
however, the normal stress in the spanwise direction was not measured. The LDV system was used to measure 
properties of the incoming boundary layers at the splitter tip. The boundary layer thickness, displacement thickness, 
and momentum thickness in each of the two streams feeding the mixing layer were measured just upstream of the 
splitter plate trailing edge. The spectra of the turbulent kinetic energy and the fundamental frequency of vortex 
shedding were not measured. 
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Top Flow: M-1.91, U=700m/s, P=49kPa 
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Bottom Flow: M=1.36, U=399m/s, P=49kPa 


Figure 1: Schematic of splitter tip separating the top flow from the bottom flow. 


Table 1. Experimental flow conditions. 


Parameter 

Top Flow 

Bottom Flow 

M 

1.91 

1.36 

U (m/s) 

700 

399 

T t ( K) 

578 

295 

P (kPa) 

49 

49 

P (kg/m 3 ) 

0.51 

0.79 


III. Formulation 

The Favre-filtered Navier-Stokes equations are solved numerically with both the high- and low-order methods. 
The viscosity and thermal conductivity are allowed to vary with temperature. Favre filtering is a density- weighted 
spatial filtering procedure that is applied to the time dependent Navier-Stokes equations to remove small-scale 
fluctuations that are too small to be resolved. The continuity, momentum, and energy equations are: 
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The density- weighted spatially-filtered velocity is denoted as u h while the spatially-filtered density and pressure are 
p and P respectively. The stress tensor is f t j. For the scales that are too small for the scheme and grid to resolve, an 
optional subgrid scale model specifies r- GS and qf GS . 


IV. Computational Approach 

Two different numerical approaches are employed herein with the purpose of investigating whether a low-order 
scheme is capable of reproducing the results of a high-order scheme. High-order methods require less grid to resolve 
turbulent structures, whereas low-order methods require additional grid points to achieve the same level of accuracy. 
However, the large stencil sizes and complex numerics typically limit high-order methods to simple geometries. If 
the low-order method is capable of producing the same answer as the high-order method, then the low-order method 
may be more attractive for other geometries that are more complex than the relatively simple shear layer that is the 
focus of this paper. 


A. Flow Solvers 
1. WRLES 

The high-order solver is the Wave-Resolving LES (WRLES) code [1-3]. It is a finite difference code which 
employs the Dispersion Relation Preserving (DRP) scheme developed by Bogey & Bailly [17-18]. The stencil size is 
11 points and is formally fourth-order accurate. While a traditional fourth-order scheme uses fewer points in the 
stencil, the additional points in the DRP scheme are used to minimize the dispersion error. It is an explicit central 
scheme similar to the classic tenth-order central scheme but with weights carefully chosen to minimize the 
dispersion error. Since it is a central differencing scheme, it is not inherently dissipative. A spatial filter, developed 
specifically to match the resolution of the differencing scheme, is used to preserve numerical stability. A coefficient 
multiplying the effect of the filter has been added, so that the magnitude of the dissipation can be controlled. The 
filter’s purpose and effects are strictly to help improve numerical stability. The DRP filter coefficient used was 0.5. 
Temporal discretization is performed using a low-dissipation and low-dispersion forth-order Runge-Kutta algorithm 
summarized in Berland et. al [18]. Three subgrid scale models have been implemented: Smagorinksy's, Vreman's, 
and a dynamic Smagorinsky. In particular, Vreman’s [22] algebraic SGS model is used herein since it is less 
dissipative than the Smagorinsky subgrid model and not as computational intensive as the dynamic model. 

WRLES uses the Message Passing Interface (MPI) standard to divide the domain across nodes and the Open 
Multi-Processing (OpenMP) directive to parallelize work on each node. In order to maintain order accuracy across 
interfaces, WRLES uses overlapping, point-matched zone coupling. 


2. Wind-US 

The low-order solver is Wind-US, the production flow solver of the NPARC Alliance. While it has historically been 
used primarily in RANS mode, the code may be used in LES or hybrid RANS-LES mode. It is capable of working 
with unstructured meshes as well as curvilinear structured meshes. It has several choices for temporal and spatial 
discretization as well as different options for flux splitting that are well documented [4] . For this particular problem, 
the numerical scheme used in the Wind-US CFD code is the second-order Roe Physical up winding finite-volume 
scheme for spatial discretization and an implicit first-order time-stepping scheme. The Roe Physical scheme is based 
upon the original Roe scheme but modified for stretched grids. Note that while implicit temporal schemes require a 
matrix inversion, they are much more stable than explicit time schemes. When run in LES mode, Wind-US employs 
the Implicit (or Numerical) LES (ILES) philosophy. As discussed by Pope [23], this approach assumes that the 
numerical scheme and grid act like a subgrid scale model, and so no explicit SGS model is used. 

For parallelization, Wind-US was run with the PVM (Parallel Virtual Machine) directive to create a parallel 
computing system. Wind-US uses abutting, point-matched zone coupling that is second-order accurate across 
interfaces. Within Wind-US, the MPI directive is also available for parallel computing. 
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B. Grid 

The baseline grid used is structured as shown in Figure 2. The baseline mesh uses 14.4 million points to resolve 
the LES mixing section and has a thin span wise width of 6 mm. Forty one points resolve the 0.5 mm splitter tip 
yielding a Ay = Ax = 0.0125 mm in this region. Grid details such as the number of points, domain size, and 
minimum and maximum grid spacing are provided in Table 2 for the mixing region. The total length of the mixing 
section is L = 500 mm as may be observed in Figure 2. Downstream of the splitter plate towards the outflow the 
grid is stretched in the x-direction. The extent of grid stretching is deduced by comparing the minimum and 
maximum grid spacing. 



nr 

48mm 

_L 


Figure 2: Boundary conditions and dimensions of the LES mixing section ( only every 8 th point is plotted). 


Table 2. Grid properties of the mixing section. 


Parameter 


Direction 



X 

y 

z 

Number of points 

1025 

425 

33 

Length of domain (mm) 

500 

48 

6 

Minimum spacing (mm) 

0.0125 

0.0125 

0.1875 

Maximum spacing (mm) 

5 

0.2 

0.1875 


C. Boundary Conditions 

1. General Boundary Treatment 

The flow conditions are taken from the experiments of Goebel & Dutton [10]. The flow in the computational 
domain is supersonic throughout, thereby simplifying inflow and outflow boundary treatments (Figure 2). Inflows 
are fixed and the outflow is extrapolated. The upper and lower tunnel walls are assumed to be sufficiently far away 
from the shear layer such that they do not affect the mixing and are set as in viscid walls. Periodic boundaries are 
imposed on the spanwise sides of the domain. Their purpose is to simulate an infinite span and avoid the cost of 
modeling the sidewalls. The experiments had a finite span of 96 mm. As mentioned previously and discussed in 
greater detail later in this paper, the width of the spanwise dimension was investigated as a significant part of this 
work. 

2. Inflow Boundary Treatment 

In the experiment, turbulent boundary layers were present on the splitter plate upstream of the mixing section. 
Dutton & Goebel [10] provided the values of the boundary layer thickness, displacement thickness, and momentum 
thickness. However, the actual velocity profiles and turbulent intensities were not reported. The best method for 
simulating these turbulent boundary layers would be a complete LES of the upstream region, producing unsteady 
turbulent boundary layers that match the experiment's integral thickness parameters. However, this type of analysis 
is cost prohibitive given the extremely small turbulent scales that would exist in a boundary layer at these Reynolds 
numbers. 

In this study, the inflow conditions were obtained from a separate two-dimensional RANS calculation of the 
mixing layer and upstream channels using the Menter Shear Stress Transport (SST) turbulence model. To do so, the 
lengths of the upstream ducts were varied until they produced RANS profiles that matched the experimental values 
of the displacement thicknesses. Table 3 compares the RANS -obtained profiles with the experiment by listing the 
computed values of displacement thickness 8*, mean streamwise velocity U , and the area-averaged pressure P. 
Notice that the RANS solution indicates that a pres sure -matched flow of 49kPa is being delivered as in the 
experiments of Goebel-Dutton [10]. The incoming velocity profiles are depicted in Figure 3 and have displacement 
thicknesses matching the experiments to within 3.5% for the top flow and within 6% for the bottom flow. The 
resulting RANS mean-flow profiles were used as fixed conditions at the inflow plane to the mixing section. A 
limitation of this boundary treatment is that the RANS method only provides a steady state profile to an unsteady 
LES mixing region [24] . 
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Table 3. Properties of the incoming flow delivered to the LES mixing section. 


Solution 

Top Flow 

Bottom Flow 


6* mm 

U, m/s 

P, kPa 

6* mm 

U, m/s 

P, kPa 

RANS Upstream 

0.9318 

700.3 

48.898 

0.4681 

398.8 

48.940 

Experiment 

0.90 

700 

49 

0.44 

399 

49 



(U-U 2 )/AU 

Figure 3: Velocity profile of the flow entering the shear layer. 

D. Time averaging and sampling procedure 

To ensure numerical stability, the maximum time-step for the implicit low-order solver was found to be dt = 
0.05 ps, while the explicit high-order solver was more stringent and the maximum time step was dt = 0.01 ps. At 
the beginning of the simulations, there is a transient wave that travels through the domain. This originates from a 
normal shock produced at solution start-up. The time for all initial transient effects to be eliminated and before 
statistics can be gathered is no less than 5.5 FTTs (Flow Thru Times). The Flow Thru Time is defined as the time it 
takes a particle inside the shear layer to travel from the splitter tip to exit: 1 FTT = L/U AVG = 0.51 milliseconds. 
U AVG is an estimate of the speed inside the shear layer taken to be the average velocity of the top and bottom 
streams. Extending the start-up wait time until 9.5 FTTs was found to be unnecessary, since statistics gathered after 
employing these two different start-up times were identical. 

The effect of sampling duration was examined by gathering second-order statistics over periods of 9.5 FTTs and 
33 FTTs that were then time averaged. Note that both of those sampling durations are considerably long. The 
resulting rms plots were compared and it was found that extending the sampling duration from 9.5FTTs to 33FTTs 
did not affect the rms values. Turbulent statistics were gathered at three frequencies corresponding to every 5 ps, 1 
ps, and 0.4 ps. Those frequencies were all extremely high and the results were found to be identical. A similar study 
varying sampling duration and frequency was performed by Mankbadi & Georgiadis [20]. After the turbulent 
statistics have been gathered, the power spectra is obtained by sampling the data every 0.05 ps after which a fast 
Fourier transform is performed to identify the fundamental vortex shedding frequency. 

V. Results and Discussion 

All the computations were performed on the NASA Pleiades supercomputing platform which is ranked amongst 
the top ten U.S. supercomputers. Pleiades is made up of tens of thousands of nodes capable of 4.49 petaflop/s. Each 
node has two processors and each processor has either 6, 8, 10, or 12 cores. A typical simulation adhered to the 
following procedure: 1) initialize the flow field and wait for the start-up transient wave to exit the domain, 2) begin 
collecting data required for computing time-averaged turbulent intensities, and 3) at a much higher sampling rate 
(about one hundred times the fundamental) gather data to compute the power spectra of the flow. This procedure 
was used for both flow solvers. Wind-US simulations typically employed 39 nodes with a total of 468 cores and the 
turnaround time for the entire procedure was about a week. WRLES utilized 129 nodes with a total of 1536 cores 
and also had a turnaround time of about a week. 

Figure 4 shows dimensionless density contours of the mixing layer obtained using the two flow solvers for the 
baseline grid with 6 mm span wise width. The high-order method, WRLES, is compared against the low-order 
method, Wind-US. Note that the size of the structures grows as the vortices pair and propagate downstream. There 
are more organized structures that can be seen in the WRLES plot whereas in the Wind-US plot the structures 
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appear to be smeared (Figure 4 A and B). The shock structure also slightly differs between the two. There are more 
reflections present in Wind-US than with WRLES. Wind-US’s small stencil and Riemann solver are adept at shock 
capturing whereas WRLES reduces the order to avoid smearing the shock with its larger stencil. The von Karman 
shedding is evident in Figure 4C and D where vortices are being produced just downstream of the splitter tip. The 
close-up views of the splitter tip demonstrate that the high-order method better resolves vortex shedding than the 
low-order method. Further downstream at x = 100 mm, the Wind-US solution appears to maintain the regular von 
Karman pattern while the WRLES solution appears to have broken down into a more irregular pattern. 
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0 5 10 x 15 20 25 0 5 10 x 15 20 25 

Figure 4: Dimensionless density contours of the mixing section for A) the high-order solver WRLES and B) low- 
order solver Wind-US. Close up views of the splitter tip showing von Karman shedding for C) WRLES and D) Wind- 


US. 
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A. Flow Solver Comparison 
1. Effects of Flow Solvers 

Streamwise mean and rms velocities from the high-order method WRLES and the low-order method Wind-US 
are compared in Figure 5 (the transverse and span wise rms velocities are presented in subsection B). Recall that the 
Wind-US approach employed a second-order finite volume solver and WRLES is a formally fourth-order accurate 
finite difference method employing an eleven-point DRP stencil. Both methods employ the ILES approach where no 
subgrid scale model is used. Despite the significant differences in numerical methodology, notice the relatively close 
agreement between the two numerical methods shown in Figure 5. The rms profiles at three different axial stations 
are very similar (Figure 5D-F), and the mean velocity profiles match fairly well (Figure 5B-C). The largest 
difference in mean velocity is in the wake at x = 50 mm in Figure 5 A. 

Figure 6 compares plots of the transverse velocity spectra of both numerical methods. The frequency is scaled 
by the fundamental frequency of vortex shedding f vs . The vortex shedding frequency, f vs , predicted is 124.51 kHz 
and 134.28 kHz by Wind-US and WRLES respectively. The data from the Wind-US solution shows a sub-harmonic 
exists at f/f vs = 0.5. The vortices produced at the tip recombine further downstream but at a lower frequency than 
that at which they were produced [25]. The vortices coalesce at the frequency of the sub-harmonic, f/f vs = 0.5. The 
sub-harmonic peak is not only present very close to the tip at x = 10mm but persists until x = 100mm. The 
fundamental frequency is the strongest at x = 10mm close to the tip and diminishes downstream, until it has 
disappeared by x = 100mm. 

The sub-harmonic cannot be seen in the spectra coming from WRLES. In addition, the spectra coming from 
WRLES shows that all the low frequency peaks have been damped out by x = 50mm and there is more energy in 
the smaller-scale, higher-frequency structures. The high-order numerics in WRLES can resolve smaller scales of 
turbulence, and these scales may be responsible for damping out the large-scale structures that persist in the Wind- 
US solution. Because of the differences in the fundamental frequency predicted and overall spectra, it is a little 
counterintuitive that two numerical methodologies would agree to such an extent in terms of the velocity and rms 
profiles shown in Figure 5. 
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Figure 5: Comparison between the low-order solver Wind-US and the high-order solver WRLES on the 14.4 million 
grid. For three different axial stations, the shear layer’s mean streamwise velocity A) at x = 50 mm, B) at x = 
100 mm, and C) at x = 150 mm. The rms streamwise velocity D) at x = 50 mm, E) at x = 100 mm, and F) at x = 
150 mm. 
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Figure 6: Spectra of the transverse velocity for both numerical methods on the 14.4 million point grid. 


2. Effects of SGS Model 

In the previous section, the effects of the numerical scheme used by the high-order and low-order methods were 
examined. This section examines the effects of subgrid scale modeling on the llpt DRP scheme solution. An 
additional case was computed with WRLES using Vreman's subgrid model. The WRLES with SGS and WRLES 
without SGS cases are compared in Figure 7 and 

Figure 8. As evident from the data, the mean profile and rms profile are unaffected. The shedding frequency 
predicted by WRLES is 134.28 kHz for both cases and only minor differences are seen in the spectra. Therefore, the 
effect of SGS modeling is negligible for the 14.4 million point grid. In Figure 9, a plot of the SGS viscosity 
normalized by the laminar viscosity is provided. It is clear that the subgrid scale model is only active within the 
shear layer. Closer inspection of the scale reveals that gi t /p does not exceed unity anywhere in the mixing section. 
This explains why SGS effects are negligible for this problem. The grid is fine enough such that the vast majority of 
turbulence is computed directly and the SGS model does not add a significant amount of dissipation. If the grid was 
coarsened, then perhaps the SGS model would come into play. The fact that p t /p does not exceed unity also helps 
explain why the second-order method performed very well. The grid was fine enough in the mixing region, making 
any benefit of a higher order scheme irrelevant. The coarser grid and SGS effects will be explored in the next 
subsection. 
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Figure 7: Effect of subgrid scale modeling on the shear layer’s streamwise mean and rms velocity at three different 
locations. The streamwise velocity A) at x = 50 mm, B) at x = 100 mm, and C) at x = 150 mm. The rms 
streamwise velocity D) at x = 50 mm, E) at x = 100 mm, and F) at x = 150 mm. 
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Figure 8: Effects of SGS modeling on the spectra of the transverse velocity. 
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Figure 9: The turbulent sub grid scale viscosity normalized by the laminar viscosity for the 14.4 million point grid. 
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3. Effect of Grid Density 

Recall in Section A-l, we demonstrated that the low-order solver could produce similar mean flow and turbulent 
statistics to the high-order solver on a fine grid. This section compares the two codes on a coarser grid. The 
coarsened grid used every other point in the streamwise and transverse directions while retaining the same spanwise 
width and points. 

The results of the grid study using the high-order solver WRLES are presented in Figure 10. Overall, the high- 
order solver performs nearly as well on the coarser grid. Close to the splitter tip, even the peak streamwise rms 
velocity values differ very little between the two resolutions. However, the mean flow predictions in the wake of the 
splitter tip are sensitive to resolution. Further away from the splitter tip, the mean flow and turbulent intensities 
predicted by the 3.6 million point grid are almost identical to the fine 14.4 million point grid. Figure 10 also 
demonstrates that effects of SGS modeling are still very minor, despite the coarser grid employed. Comparing the 
SGS viscosities of Figure 9 for the finer grid to those of the coarsened grid in Figure 11, demonstrates that as the 
resolution was coarsened the SGS viscosity increases, but is still on the same order as the laminar viscosity. 

The results of the grid study using the low-order solver Wind-US are presented in Figure 12. These results 
demonstrate that the low-order method changes significantly at the coarse resolution. The low-order method cannot 
resolve the complex flow physics of vortex shedding, mixing, and shear layer growth. The inherent dissipation of 
the low-order scheme dampens out the instabilities and inhibits proper turbulent growth. There are substantial 
differences in both the mean flow velocities and the rms quantities. Not only is the peak rms value diminished but 
the extent of the profile across the shear layer is reduced significantly from that observed with the higher resolution. 
In other words, very little turbulent mixing is occurring. 

Figure 13 shows the transverse velocity spectra comparing the high-order and low-order methods at the coarse 
resolution. Note that the high-order method has well-defined sharp peaks at the fundamental frequency and its 
multiples. However, the low-order method has a blunt peak at the fundamental frequency. Recall the fundamental 
frequency predicted by WRFES on the fine grid was 134.28 kHz, here for the coarse grid it is higher at 141.91. 
Fikewise recall Wind-US predicted a fundamental frequency of 124.51 kHz, for the coarser grid it is also higher 
128.17 kHz. 

These findings demonstrate the expected result that high-order algorithms need less resolution than low-order 
methods to achieve the same level of accuracy. Unlike high-order methods, low-order methods are expected to fail 
at coarse resolutions and indeed did in this study. The increase in accuracy that high-order methods offer means that 
high Reynolds number problems become tractable at more reasonable grid sizes. With regards to computational wait 
time, high-order methods typically incorporate multi-stage time stepping and use large stencils that drives up the 
computational cost. Nevertheless, in many cases high-order methods can deliver that same level of accuracy at a 
cheaper computational cost since less grid is required. However, for the current problem, the two methods’ 
computational costs were about the same. This is primarily due to two reasons: (1) the explicit time-stepping of the 
high-order code requires a smaller time step than the implicit time stepping of the low-order method, and (2) the 
low-order method needs to regularly issue a checkpoint where the entire data is written to disk then an external 
process strips down the file to keep only the desired data. If not for this limitation in I/O performance, the low-order 
method would likely have been much faster than the higher-order method for the same grid. 
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Figure 1 0: Grid study using WRLES with a spanwise width of 6mm. The mean streamwise velocity A) at x 
50 mm, B) at x = 100 mm, and C) at x = 150 mm. The rms streamwise velocity D) at x = 50 mm, E) at x 
100 mm, and F) at x = 150 mm. 


15 



0 50 100 150 

X 

Figure 11: The turbulent subgrid scale viscosity normalized by the laminar viscosity for the 3.6 million point grid. 
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Figure 12: Grid study using Wind-US with a spanwise width of 6mm. The mean streamwise velocity A) at x 
50 mm, B) at x = 100 mm, and C) at x = 150 mm. The rms streamwise velocity D) at x = 50 mm, E) at x 
100 mm, and F) at x = 150 mm. 
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Figure 13: Spectra of the transverse velocity for both numerical methods on the 3.6 million point grid. 
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B. Effect of Spanwise Width 


Table 4. Spanwise width study using WRLES with SGS. 


Spanwise Width 

Resolution 

Grid Points (million) 

Max Wrms /AU 


6mm 

1025 x 425 x 33 

14.4 

x = 50 mm 
0.0947 

100mm 

0.09343 

150mm 

0.07341 

12mm 

1025 x 425 x 65 

28.3 

0.1005 

0.1248 

0.1139 

24mm 

1025 x 425 x 129 

56.2 

0.1056 

0.1362 

0.1315 

48mm 

1025 x 425 x 257 

112.0 

0.1106 

0.1356 

0.1253 


Periodic boundary conditions are imposed in the spanwise direction to simulate an infinite span and avoid 
modeling the sidewalls. Selection of too small of a spanwise width cannot guarantee accurate capturing of the three- 
dimensional nature of vortex dynamics. To avoid imposing a nonphysical spanwise mode, the spanwise extent and 
resolution must be large enough to capture the wavelengths of the largest turbulent structures in the domain. 

An insufficient spanwise depth has been shown to artificially dampen the spanwise formation of vortices and 
resemble a 2D instability [20]. In some cases, a 2D instability predicts higher growth rates than a three dimensional 
instability [26] . It is necessary to have a sufficient spanwise width to allow for the formation of three-dimensional 
coherent structures. Four simulations are performed here using the high-order method, with various spanwise widths 
in order to determine if: (1) artificial suppression is present when the spanwise width is smaller than the length scale 
of the largest coherent structure, (2) this artificial suppression is eliminated as the spanwise width is extended to be 
several times larger than the dominant length scale, and (3) as the mixing layer grows downstream, the size of the 
coherent structure grows thereby requiring a larger spanwise depth to be resolved. 

WRLES with the SGS model was used to study the effect of the spanwise width. Table 4 shows the grid 
resolution as the spanwise domain is doubled three times consecutively. Georgiadis et al. [19] simulated this same 
flow and noted significant differences in rms statistics when examining much smaller spanwise depths ranging from 
lmm to 6mm. The baseline grid was chosen to have the largest spanwise depth considered by Georgiadis et al. [19] 
of 6 mm and had a spanwise grid spacing of A Z = 6mm/32. Keeping the spanwise grid spacing fixed, the spanwise 
width was extended up to 48mm, utilizing a total of 112 million grid points. The peak values of W rms were also 
obtained. The streamwise and transverse resolution was unchanged. 

Figure 14 demonstrates the effects of increasing the spanwise width by plotting the streamwise vorticity at the 
three axial stations for spanwise widths of 6 mm, 12mm, 24mm, and 48mm. Evidently the size of the structures is 
small at x = 50mm, but as they propagate downstream the structures grow in size and by x = 150mm the 
structures have grown considerably. The effect of increasing the spanwise width is most easily observed at x = 
150mm. For a spanwise width of 6mm, the development of the structures is suppressed, but as the spanwise width 
is increased to 48mm, the structures are allowed to develop and thus appear to be bigger. The size of the turbulent 
structures relative to the mixing layer thickness can be seen by comparison with the value of b that is shown to the 
right of Figure 14. 
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Figure 14: Effect of spanwise width on the dimensionless streamwise vorticity 
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Although the mean flow is not affected by spanwise widths, the turbulence intensities are sensitive. The 
turbulence intensities in all three directions and turbulent kinetic energy are compared for all four spanwise widths at 
x = 50 mm (Figure 15), x = 100mm (Figure 16), and x = 150mm (Figure 17). At axial stations of 50mm and 
100mm, the peak values of u rms and v rms are almost unchanged as the spanwise width was gradually increased 
from 6mm to 48mm. At x = 150mm (Figure 17A-B), the peak value in v rms is also unchanged but the peak u rms 
value slightly increases. From examining u rms and v rms , the need to extend the spanwise domain is not apparent 
since u rms and v rms are not changing appreciably. However, close inspection of w rms and k shows that w rms is 
artificially being suppressed for the 6mm case (Figure 15C-D, Figure 16C-D, Figure 17C-D). Expanding the width 
of the domain increases the peak w rms at the x = 50mm and x = 100mm locations and k correspondingly 
increases. This differs from a similar study of spanwise width for a square cylinder computation [20]. For that 
problem w rms also increased with an increase in spanwise width but the k remained constant, forcing a rebalance of 
the energy between u rms> v rms , and w rms . 

Note that at x = 50mm, where the vortex size is small, the effect of the width on w rms is relatively small. 
However, as the turbulent structures increase in size moving further downstream, a larger spanwise width is 
required. Hence, at x = 100mm, note that the peak rms value is identical between the 24mm and 48mm case. 
Further downstream at x = 150mm, the peak value of w rms of the 48mm case lies between the 12mm and 24mm 
cases. This demonstrates that extending the spanwise width eliminates artificial suppression of the spanwise rms 
value. 

In the aforecited work of Mankbadi & Georgiadis [20] on a square cylinder, a spanwise depth ten times greater 
than the cylinder’s diameter demonstrated the W rms convergence. For a cylinder at low Reynolds number, the shed 
vortices are about the size of the cylinder’s diameter. However, for the shear layer examined here, the relevant 
length scale is the shear layer width, which grows downstream. The Kelvin-Helmholtz instability sets off the vortex 
shedding process that then becomes turbulent mixing downstream. As a consequence, the size of the largest 
structure increases axially. Hence, at x = 50mm where the length scale is small, a spanwise depth of 6mm is 
sufficient and spanwise width convergence is observed in Figure 15C. Further downstream at x = 150mm, the size 
of the length scale has grown and hence a variation in w rms is observed in Figure 17C. 

The relevant length scale is the shear layer thickness which correlates with the size of the largest coherent 
structure. A closer examination of the shear layer thickness is necessary to assess the ratio between the thickness and 
spanwise width at which spanwise width independence is reached. The shear layer thickness was defined by Goebel 
& Dutton [10] as b = y 0 . 90 — yo.io> where y 0 90 is the location corresponding to U* = 0.90. Recall U* is normalized 
such that at the top flow it is one and at the bottom flow it is zero. From Table 5, note that at x = 50mm the shear 
layer thickness is b = 2.5mm, while it has grown to the shear layer thickness of 4.4mm at x = 150mm (see Table 
5). 

Note that the U rms for the 6mm case in Figure 15 A is identical to the 48mm case, thus spanwise width 
independence has been reached at 6mm. When the 6mm width is normalized by the 2.5mm thickness at x = 
50mm it yields a ratio of 2.3. Likewise in Figure 17C the peak w rms value for the 48mm case lies between the 
12mm and the 24mm cases, suggesting that they are relatively close at a spanwise width of 12mm. When the 12mm 
width is normalized by the 4.4mm thickness at x = 150mm it yields a ratio of 2.7. These results suggest that the 
solution becomes domain independent when the spanwise width is 2.5 times larger than the mixing layer thickness. 

For the shear layer case investigated herein, it was possible to extend the domain out until 48mm. This should 
ensure uninhibited growth of the turbulent structures, even at the furthest downstream measurement location of 
450mm. For stations much beyond that, the spanwise width would likely need to be extended even more. However, 
at some point the spanwise domain cannot be extended without taking sidewall effects into consideration. The 
proposed factor of 2.5 cannot necessarily be generalized to other more complex flows. For other similar problems 
where resolution is scarce, a spanwise width that is only 2.5 times larger than the relevant length scale may be 
sufficient to prevent artificial suppression of W rms . 
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Figure 15: At x = 50 mm, the effect of spanwise width on the shear layer’s A) streamwise rms velocity, B) 
transverse rms velocity, C) spanwise rms velocity, and D) the turbulent kinetic energy. 
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Figure 16: At x = 100 mm, the effect of spanwise width on the shear layer’s A) streamwise rms velocity , B) 
transverse rms velocity, C) spanwise rms velocity, and D) the turbulent kinetic energy. 


22 


0.4 

0.3 

0.2 

0.1 

t o 
- 0.1 
- 0.2 
-0.3 
-0.4 


|x=150mi 

m (A) 

Y • 






/ f / • 

/ • 

/// • 



6mm 

■ 12mm 
24mm - 

!// • 



48mm 




• EXP 


0.05 


0.1 0.15 

u /All 


0.2 


0.25 


0.4 

0.3 

0.2 

0.1 

0 

- 0.1 

- 0.2 

-0.3 

-0.4 


0.05 


x=150mi 

m 

(C) 













— 6mm ■ 
12mm 




. 24mm 




48mm 


0.1 0.15 

W /All 


0.2 


0.25 




k/AU 2 


Figure 17: At x = 150 mm, the effect of spanwise width on the shear layer’s A) streamwise rms velocity, B) 
transverse rms velocity, C) spanwise rms velocity, and D) the turbulent kinetic energy. 


Table 5. Shear layer thickness. 


Location 

LES (z = 48mm) 

Exp. 

x (mm) 

b(mm) 

b(mm) 

50 

2.549 

3.241 

100 

3.480 

4.508 

150 

4.426 

5.284 
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VI. Conclusions 

This work compared a high-order method with a low-order method for performing large-eddy simulations 
(LES). The high-order method (WRLES) employed an eleven-point DRP scheme that is formally fourth-order 
accurate in space. For time stepping, it used the four-stage Runge-Kutta method, which is fifth-order accurate 
locally. For subgrid scale modeling, it uses Vreman’s model. The low-order method (Wind-US) is second-order 
accurate spatially and uses the Roe Physical scheme. For time stepping, it used an implicit backward Euler step that 
is first-order accurate locally, and does not utilize any subgrid scale modeling. EES were performed using both 
methods to capture the turbulent mixing between two supersonic streams at a convective Mach number of 0.46. The 
high-order and low-order methods were evaluated at two different levels of grid resolution. In order to compare with 
the experiment, the inflow profiles for LES of the shear layer were obtained using Reynolds Averaged Navier- 
Stokes (RANS) modeling over the upstream splitter plate such that the displacement thickness matched the 
experimental value reported. The results presented in this paper are summarized as follows: 

• For a fine grid resolution, the low-order method produced a similar solution to the high-order method. Both 
methods predict turbulent stresses that are in good agreement with experimental data. 

• Effects of subgrid scale modeling were found to be negligible at both resolutions using the higher order 
method and the Vreman SGS model. At even coarser resolutions, the SGS model might have a larger 
effect. 

• WRLES and WIND-US density contours were compared. For both methods, vortex shedding was observed 
to originate at the splitter tip and quickly transition to turbulence. More organized structures were seen with 
the high-order method whereas those structures appeared to be smeared with the low-order method. 

• The time for all initial transient effects to be eliminated and before statistics can be gathered is no less than 

5.5 FTTs (Flow Through Times). The recommended sampling duration is 9.5 FTTs or longer. Any time- 
averaging window that meets these two conditions will not produce any substantial uncertainty. 

• When the DRP filter was nearly turned off, it was affirmed that its purpose is primarily for numerical 
stability and does not act as a SGS model. Thereby confirming that the reason the low-order and high-order 
method agreed was because the grid was sufficiently fine enough. 

• For a coarse grid of 3.6 million points, the high-order flow solver arrived at the same solution it reached 
with the 14.4 million point fine grid. However, the low-order solver did not demonstrate grid independence 
between the coarse and fine grid. Thus, high-order algorithms need less resolution than low-order methods 
to achieve the same level of accuracy. 

• We conducted a parametric study for the effect of the span wise width on the accuracy of the solution by 
varying the width from a 6mm depth to a 48mm depth. As the shed vortices travel downstream and grow in 
size, the spanwise width may not be large enough to resolve all the spanwise modes and may introduce 
artificial ones. The peak root-mean- square of the spanwise fluctuation, w rms , was found to increase by 
fifty percent, when the spanwise width was increased from 6mm to 48mm. A spanwise width that is at least 

2.5 times larger than mixing layer thickness was found to be necessary to capture three-dimensionality of 
vortex shedding. 

• Despite the differences in parallelization and time-stepping scheme, the two method’s computational costs 
were about the same. The larger stencil size and explicit time-stepping of the high-order code resulted in a 
larger computational cost for a given grid size compared to the low-order, implicit time-stepping code. 
However, the input/output performance of the low-order code was slowed by the need to regularly 
checkpoint the solution and spawn an external tool to extract the unsteady flow information. 

The high-order method herein proved to be more accurate on a coarse mesh when compared against the low- 
order method. The small stencils of low-order methods are advantageous for their ability to handle complicated 
geometries whereas high-order methods are constrained by a large stencil size. The work herein establishes that low- 
order methods can be used for LES but require a larger number of grid points than their counterpart. Historically, 
most high-order methods do not retain their accuracy with unstructured grids. Recently, unstructured grids are 
favored for the inherent ability to redistribute mesh points where needed and remove unnecessary mesh points. 
Ideally the benefits of low-order and high-order methods would be combined by creating stable high-order 
unstructured algorithms. Such a method may push the state of the art by enabling high-order unstructured LES of 
practical applications. 
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